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ABSTRACT 



Context. In isothermal disks the migration of protoplanets is directed inward. For small planetary masses the standard type I migration 
rates are so fast that this may result in an unrealistic loss of planets into the stars. 

Aims. We investigate the planet-disk interaction in non-isothermal disks and analyze the magnitude and direction of migration for an 
extended range of planet masses. 

Methods. We have performed detailed two-dimensional numerical simulations of embedded planets including heating/cooling effects 
as well as radiative diffusion for realistic opacities. 

Results. In radiative disks, small planets with M p i anel < 50M Earth do migrate outward with a rate comparable to absolute magnitude of 
standard type I migration. For larger masses the migration is inward and approaches the isothermal, type II migration rate. 
Conclusions. Our findings are particularly important for the first growth phase of planets and ease the problem of too rapid inward 
type-I migration. 
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1. Introduction 

Planets form in disks surrounding young stars. The growing 
protoplanets undergo an embedded phase where the gravita- 
tional interaction with the ambient gaseous disk results in a 
change of its orbital elements. For protoplanets with masses be- 
low about 30 Earth masses the disk is not disturbed too strongly 
and the interaction can be treated in the linear approximation. 
Calculations of the total disk torques acting on the planet lead 
generally for these small masses to a reduction of the semi- 
major axis, i.e. to an inward migration dGo ldreich & Tremamel 
[1971 IWardl Il997t lTanaka"etllTl200l iTanaka & Wardl 12001 . 
It soon turned out that the inward drift of this type-I migra- 
tion is very fast and the planets might be lost befor e they 
can grow to larger objects (Korycansky & Pollack 1993). This 
problem has become more visible after comparing population 
synthesis mod els with the chara c teristics of obser ved plane- 
tary systems dAlibert et al l l2004t llda & Linl 120081) . To avoid 
this rapid phase of inward migration, alternative scenarios have 
been sought. In a turbulent disk, migration occurs stochasti- 
cally with inward and o utward phases which slows down the 
migration dNelsonl 12005). Departures from the linear regime at 
around 10-20 Mnarth can also lead to reduced inward migration 
(Masset et al. 2006a). However, both processes are not sufficient 
to solve th e problem. The pla net trap scenario to halt planetary 
migration (Ma sset et al.ll2006bT) requires a positive density gra- 
dient which may not be given in general. 

To simplify the calculations, nearly all of the analytical 
and numerical studies devoted to study the planet-disk in- 
teraction process have focussed on isothermal disks, where 
the temperature is a given function of the position in the 
disk. Early work on non-isothermal disks focussed on high 
mass embe dded planets and did not notice a strong e ffect on 
migration (D'Ang elo etafl 120031: iKlahr & Klevl 120061) . Using 
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a fully three-dimensional r adiative calculations of a n em- 
bedded small mass planet, Paardekooper & Mellemal |2006) 
have shown in a very important work that migration can 
be significantly slowed down or even reversed when ther- 
mal effects are included. Subsequent analysis indicate that 
this beha viour is related to a radial entropy gradient in 
the flow (Baruteau & Masset 2008; Paardekooper & Mellema 
2008). Recentlv. lPaardekooper & Papaloizoul ( l200 8) have shown 
that for small planet masses, the combination of radiative and 
viscous diffusion may allow for long-term unsaturated positive 
torques and possible outward migration. 

In this letter we investigate this possibility in more detail for 
a whole range of planetary masses, which will allow us to es- 
timate its effect on the long term evolution of the planet. For 
that purpose we perform two-dimensional numerical hydrody- 
namical simulations of embedded planets in radiative disks. A 
method to treat the three-dimensional radiative ttansfer approx- 
imately in these 2D simulations will be outlined in the next sec- 
tion. Our results on the migration rate for various masses (in 
Section 3) indicate that for masses smaller than about 50 AfEarth 
the torques remain unsaturated in the long run and migration is 
indeed directed outward, while larger planets drift inward. The 
consequence for the migration process and the overall evolution 
of planets in disks is discussed. 



2. Physical modelling 

The protoplanetary disk is treated as a two-dimensional, non- 
self-gravitating gas that can be described by the Navier-Stokes 
equations. The embedded planet is modelled as a point mass that 
orbits the central star on fixed, circular orbit. For the planetary 
potential we use a smoothing length of e = 0.6H where H is the 
vertical scale height of the disk. To calculate the gravitational 
torques acting on the planet we apply a tapering function to ex- 
clude the inner parts of the Hill s phere of the plane t, where the 
transition lies at 0.8 Hill radii (see lCrida et~aT1l2008l Fig. 2). 
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2.1. Energy Equation 

In the present setup we utilise fully radiative models with an im- 
proved thermodynamic treatment using the thermal energy equa- 
tion in the following form 



<9£c v r 



dt 



+ V • (Zc v Tu) = -pV u + D-Q- 2HV ■ F. 



(1) 



Here u = (u r , u^,) is the two-dimensional velocity, 2 the den- 
sity, p the pressure, T the (mid-plane) temperature of the disk, 
and c v is the specific heat at constant volume. On the right hand 
side, the first term describes compressional heating, D the (ver- 
tically averaged) dissipation function, Q the local radiative cool- 
ing from the (vertical) surfaces of the disk, and F denotes the 
two-dimensional radiative flux in the (r, ^)-plane. For all mod- 
els H is calculated from the sound-speed as H(r) = c s / Sl^ir), 
where Qk is the Keplerian angular velocity, and c s = yJypfL, 
y — 1 .43 being the adiabatic index. The 2D-pressure is given by 
P = Rgos^T Ip, with the mean molecular weight p =2.35. 

To calculate the radiative losses Q (fro m the two sides o f the 
disk) we follow lD'Angelo et all (120031) and lKlev et ail d2005l) . 



Q = 2cr R r, 



eff ' 



where <x R is the Stefan-Boltzman n constant and T e g is an esti- 
mate for the effective temperature (lHubenvllT9 90). given by : 



with 



V3 J_ 
4 + 4t 



For our two-dimensional disk we approximate the mean vertical 
optical depth by r = (1/2)kE, wher e for the Rosseland mean 
opacity k the analytical formulae by Lin & Papaloizou (1985) 
are adopted. The radiative transport in the plane of the disk is 
treated in the flux-limited diffusion approximation where the flux 
is given by : 



Ac4aT 3 

PK 



■ vr. 



Here c is the speed of light, a the radiation const ant, p - Y,/ (2H) 
the mid-plane density, and A the flux-limiter (see Kley 1989))- 

In the following we present results of numerical simulations 
using various sub-parts of the energy equation, Eq. ©. If only 
the first term on the right hand side is used, the model is called 
adiabatic ; when only the last term is omitted, it is a model with 
local heating and cooling ; and the usage of the full energy equa- 
tion is termed fully radiative. To make contact with previous re- 
sults we will also use isothermal models where a pre-described 
radial temperature distribution is held fixed and no energy equa- 
tion is solved at all. Please note that in the full version, the ra- 
diative treatment (in Eq. ([l} ) incorporates full 3D effects of the 
radiative transport in that the vertical part (z-direction) is taken 
care of by the local cooling term, Q, and the horizontal part 
through F. 

2.2. Reference Model 

The two-dimensional (r - <p) computational domain consists of 
a complete ring of the protoplanetary disk centred on the star, 
extending from r m j n = 0.4 to r max = 2.5 in units of ro = a j up - 
5.2AU. The mass of the central star is one solar mass, and the 
total disk mass inside [r m i n , r max ] is Mdisk - O.O1M . For the 
present study we have kept the kinematic viscosity constant with 
v = 10 I5 cm 2 /s, a value which relates to an equivalent a at ro of 
0.004 for a disk aspect ratio of H/r = 0.05, where v = aH 2 Q.K- 



To construct a joint reference model for the different types 
of approximations to the energy equation (isothermal, adiabatic, 
etc.) we construct a fully radiative model with no embedded 
planet, using the above layout. 

This model is obtained from an approximate initial state 
(with H/r = 0.05) that is then relaxed to its equilibrium by 
performing a time evolution of the disk solving the full Navier- 
Stokes equations including the energy. The initial disk stratifica- 
tion at the start of this process is given by E(r) oc r~ l/2 , T(r) oc 
r _1 , and pure Keplerian rotation (u r — 0, u v = (GM t /r) 1 ^ 2 ). The 
evolution towards the equilibrium takes about 150 orbits and 
the resulting density and temperature distribution is displayed 
in Fig.Q] The density has a E oc r -1 ' 2 profile which follows for 
constant v directly from the angular momentum equation. The 
temperature has a T oc r 1 6 profile which follows from Q — D 
for the used opacity. These relations are superimposed in Fig.Q] 
For small radii (higher temperatures) the opacity law has dif- 
ferent temperature dependence and the slope T(r) changes. The 
relative scale height for this radiative model is not constant but 
falls off with radius. At r — 1 we find H/r as 0.045 and at the 
outer radius H/r « 0.03. 
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Fig. 1. Density and temperature vs. radius for the relaxed refer- 
ence model. For illustrative purpose, simple comparison power 
laws have been added. 



2.3. Numerics 

We work in a rotating reference system, rotating with the or- 
bital frequency of the planet. For our standard cases we use an 
equidistant grid in r and <p with a resolution of 128 x 384. In an 
effort to ensure a uniform environment for all models and mini- 
mize disturbances (wave reflections) from the radial boundaries 
we impose at r^ and r max damping boundary conditions where 
both velocity components are relaxed towards their initial state 
on a timescale of the order of the local orbital period. As the 
initial radial velocity is vanishing, this damping routine ensures 
that no mass flows through the boundaries such that the total disk 
mass remains constant. The angular velocity is relaxed towards 
the Keplerian values. For the density and temperature we apply 
closed radial boundary conditions. In the azimuthal direction, 
periodic boundary conditions are imposed for all variables. 

The numerical details of the used finite volume code (RH2D) 
releva nt for these disk simulations have been described in iKlevI 
(1 19991) . where we h ave additionally implemented the FARGO 
method dMassetl2000l) that speeds up the numerical computation 
of differentially rotating flows. The energy equation Eq. ([TJ is 
solved explicit-implicitly applying operator-splitting. The heat- 
ing and coolin g term D — Q 'is tre ated as one sub-step in this 

~iDEP 



procedure (see Giinther et al.ll2004l) . and the additional radiative 
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Fig. 2. Time evolution of the specific total torque exerted by the 
disk on an embedded protoplanet of 20ME ar th for various approx- 
imations to the energy equation. 



diffusion part in the energy equation is solved through an im- 
plicit method to avoid possible time step limitations. 



3. Results 

3.1. A model with 20 M Earth 

To illustrate the influence of the 4 different formulations of the 
energy equation on the torque, we present a model with an em- 
bedded 20 ME art h planet. For this particular mass, a recent 3D- 
study to analyze the effects of an isothermal disk on all orbital 
elements, indicated negative rates f or migration as well a s for 
eccentricity and inclination changes (ICresswell et al.l l2007). The 
starting configuration (at t — 0) for the 4 cases has been the 
equilibrated reference model corresponding to the fully radia- 
tive case as described above (see Fig.[T). In Fig.[2]we display the 
time evolution of the torque acting on the planet, for the different 
formulations of the energy equation. As expected, the isothermal 
condition leads after about 50 orbital periods to a constant nega- 
tive torque implying inward migration. The adiabatic model has 
an initial phase of positive torques, before settling to negative 
values of about 40 % of the isothermal case. This behaviour has 
been suggested by Baruteau & Masset (2008) who argue that in 
the adiabatic case the torques will be unsaturated and positive 
during the onset phase, while long-term calculations should con- 
verge to negative values. 

In contrast, both radiative models settle to constant positive 
values. Here, the fully radiative case yields a 25% smaller value 
due to the included heat diffusion in the disk plane. In all cases, 
we have continued the models with double and quadruple resolu- 
tion. In general, we find (for this planet mass) similar results and 
a convergence of the torques at about double resolution. Only the 
model with pure heating/cooling (i.e. no radiative diffusion) ap- 
parently requires much higher resolution for convergence, a fact 
we attribute to the highly local character of this type of energy 
equation. The addition of radiative diffusion in the disk plane 
eases the numerical requirements and makes at the same time 
the simulations physically more realistic. 

3.2. Variation of the planet mass 

To estimate the influence of the planet mass we performed 
a sequence of models with masses varying from q = 10~ 5 
(3MEanh) up to q = 10~ 3 (Mjupite,). The resulting equilibrium 
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Fig. 3. Comparison of the specific torque exerted on an em- 
bedded protoplanet for the isothermal and fully radiative model 
(thick lines). The lower panel is an enlargement of the upper. 
Two analytic curves are superimposed (thin lines). 



specific torques are displayed in Fig.[3]for the fully radiative and 
isothermal case. Clearly, for small planetary masses up to about 
50MEanh the torque exerted on the planet is positive for the radia- 
tive model and turns negative above this mass value (see bottom 
panel). On the contrary, for the isothermal models the torques are 
negative throughout and follow the expected well known results. 
From the superimposed dotted line, it is clear that the specific 
torque is proportio nal to q, as predicted by linear theory (eg. 
iTanaka et al.ll2002l) . 

For large planets (q > 2.5 x 10~ 4 ), there is little difference 
between radiative and isothermal di sks, and the torque is well ap - 
proximated by the model given bv lCrida & Morbidellil d2007l) : 
the dot-dashed curve in the top panel comes from their Eq. (15), 
with v r replaced by | * This model refines the type II migration 
rate for planets whose gap is not completely empty ; therefore, 
it applies here for q > 3 x 10~ 4 (see the gap profiles in Fig . [5j ■ 
We find that this model, and type II migration in general, is also 
valid in radiative disks. 

The two-dimensional density distribution for the two cases 
(isothermal and fully radiative) is displayed in Fig.@]for a plane- 
tary mass of q = 10~ 4 (33M £flrt /,) which has the maximum posi- 
tive torque. In the radiative model the higher temperature results 
in slightly larger opening angles of the spiral arms with a re- 
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Fig. 4. Two-dimensional grey-scale plot of the density (oc En- 
sealing between 150 and 350 g/cm 2 ) in the vicinity of the 
planet for the isothermal model (left) and the fully radiative 
model (right). The planet mass in this case refers to 1O~ 4 M or 
33MEarth- For these plots we have used double resolution models 
(256 x 768). 



duced density, while the isothermal situation allows for a higher 
mass concentration within the Roc he-lobe of the planet. As 
expected by previous investigatio ns dBaruteau & Massetl 120081: 
Paardekooper & Papaloizou 2008), the inner half of the horse- 
shoe region ahead of the planet (in Fig. [4] below the planet) is 
denser in the radiative case than in the isothermal one, due to 
the incoming of colder gas from the outer half after an U-turn. 
Symmetrically, the outer half of the horseshoe region is slightly 
depleted behind the planet. As a consequence of these changes in 
density, the Lindblad torques are slightly reduced and the coro- 
tation torque is more positive. As the net result we find for this 
case an outward migration. To sustain the necessary unsaturated 
corotation torque on the long run, the libration timescales have 
to be comparable to the cooling tim es dBaruteau & Massetl2008t 
iPaardekooper & Papaloizou 2008). For our parameters we find a 
cooling time t coo ; = TYc v jQ of about 200 Q^ 1 at the location 
r p of the planet. This is comparable to the horseshoe libration 
time, Tut, = 8^t p /(30jcjc s ) with the half-width of horseshoe- 
orbit, x s = 1.16/> ^q/(H/r) dBaruteau & Massetl [2008b . In our 
case, for q = \0 A we find Tig, = 1500^'. 




Fig. 5. Azimuthally averaged density for isothermal and radia- 
tive cases for different planet masses. 

A comparison of the radial density stratification of the 
isothermal and radiative cases (Fig. |5]l indicates shallower gaps 
for the latter, in particular the region inside of the planet (r < 1) 
is less cleared. Larger planet masses lead to gap opening with 
similar gap depths for the isothermal and radiative case. For the 
largest planet mass the gap seems to be wider in the isother- 
mal case, possibly due to lower temperatures, but boundary ef- 



fects (at r,„i„ and r max ) begin to become visible. For larger planet 
masses a radially increased domain is clearly required. 



4. Summary 

We have performed an investigation of the migration of planets 
in disks using two-dimensional numerical simulation including 
heating/cooling effects as well as radiative diffusion. 

Using different formulations of the energy equation, we 
first show that for a planet mass of 20ME art h, migration is di- 
rected inwards in the isothermal and adiabatic situation, while 
inclusion of radiative effects leads to an outward migration of 
the planet. This finding supports the torque reversal mecha- 
nis m in radiative disks due t o corotation effects as suggested 
by iBaruteau & Massetl d2008l) . A detailed parameter study for 
planetary masses in the range between 10~ 5 and 10~ 3 M o shows 
that the effect is limited to planets in the low mass regime, 
M p < 50ME a rth where corotation effects are indeed important. 
Larger mass planets open up gaps in the disk and the migration 
rate becomes similar to the isothermal case. 

Our findings are particularly important for the first growth 
phase of planets and may ease the problem to too rapid inward 
type I migration. Depending on the mass accretion rate onto the 
planet the growing planetary embryos can spend an extended 
time span in an outward migration phase and avoid loss into 
the star. However, close-in planets exist for a range of plane- 
tary masses, according to the observations. Thus, a significant 
and long outward migration phase may create new difficulties. 
Whether this problem really exists can only be answered by fol- 
lowing the actual long-term migration of planets through the disk 
including its mass growth. 

Constructing the necessary migration histories of planetary 
cores, to be used in population synthesis models, requires suit- 
able scaling laws for the migration process as a function of disk 
parameter (X(r), T{r)) for realistic accretion disks with net mass 
flow. The present study can be used as a starting point for these 
larger parameter studies. The inclusion of three-dimensional ef- 
fects and additional physics (MHD, self-gravity, mass accretion) 
will make the models even more realistic in the future. 
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